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Abstract. Certain dust particles in space are expected to appear as clusters of individual grains. The morphology 
of these clusters could be fractal or compact. In this paper we study the extinction by compact and fractal 
polycrystalline graphitic clusters consisting of touching identical spheres, based on the dielectric function of 



graphite from Draine & Lee (1984). We compare three general methods for computing the extinction of the 



clusters in the wavelength range 0.1 — 100 /im, namely, a rigorous solut ion (G erardy & Ausloos 1982) and two 



differ ent discrete-dipole approximation methods - MarCODES (Markel 1998) and DDSCAT (Draine & Flatau 



1994| ). We consider clusters of = 4, 7, 8, 27, 32, 49, 108 and 343 particles of radii either 10 nm or 50 nm, 
arranged in three different geometries: open fractal (dimension D = 1.77), simple cubic and face-centred cubic. 
The rigorous solution shows that the extinction of the fractal clusters, with A < 50 and particle radii 10 nm, 
displays a peak within 2% of the location of the observed interstellar extinction peak at ~ 4.6/im~^; the smaller 
the cluster, the closer its peak gets to this value. By contrast, the peak in the extinction of the more compact 
clusters lie more than 4% from 4.6f^m~^ . At short wavelengths (0.1 — 0.5/xm), all the methods show that fractal 
clusters have markedly different extinction from those of non-fractal clusters. At wavelengths > 5^m, the rigorous 
solution indicates that the extinction from fractal and compact clusters are of the same order of magnitude. 
It was only possible to compute fully converged results of the rigorous solution for the smaller clusters, due to 
computational limitations, however, we find that both discrete-dipole approximation methods overestimate the 
computed extinction of the smaller fractal clusters. 

Key words. Methods: numerical - scattering - dust, extinction - ISM: general 



1. Introduction 

The shape of interstellar and circumstellar grains is still 
an outstanding issue. For many years, the complexity of 
the electromagnetic scattering problem to solve has lim- 
ited the shapes studied to spheres, infinite cylinders and 
spheroids. However, the shape of many interstellar grains 
are expected to be non-spherical and maybe even highly 
irregular. One way to deal with irregular particles and so 
with clusters of dust grains is to assume that they con- 
sist of touching spheres. With such an assumption it is 



possible to construct many distinct different morphologies 
which can then be compared with observations. 

The problem of evaluating the extinction efficiency 
(Qcxt) is that of solving Maxwell's equations with appro- 
priate boundary conditions at the cluster surface. A so- 
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lution was formulated by Lorenz (189C) and Mie (1908) 
for a homogeneous single sphere and the complete formal- 
ism is therefore often referred to as Lorenz-Mie theory. 
A complementary solution bas ed on expansion of scalar 
potentials was given by Debye ( 1909 ). A detailed descrip- 
tion of this exact electromagnetic solution can be found 
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in Bohren & Huffman (1983). For a comprehensive review 
on the optics of cosmic dust see Voshchinnikov ( |2002| ). 

An updated overview of exact theories and numeri- 
cal techniques for computing the scattered electromag- 
netic field by clusters of particles is given in (Mishchenko 
et al. pOOOal ; |2000bt Fulle r & Mackowski |2000t Ciric & 
Coorey 200C; Draine [2000 ; Voshchinnikov 2002 and refer- 
ences therein). All of these methods are based on solving 
Maxwell's equations. For clusters of spheres embedded in 
non-absorbing media one such method is based on the 
Gerardy and Au sloos theory (Gerardy & Ausloos 198C; 
1982; 1983| ; 1984) which recently has been extended to also 



treat cl usters embedded in absorbing media (Lebedev & 
Stenzel |l999t Lebedev et al. |l999|) . Another method that 
is often used in practice is the discrete dipole approxima- 
tion (DDA); it has been used in a wide range of scattering 
problems concerning clusters of particles including the ex- 



tinction of inter stella r dust grains (e.g. Vaidya et al. 2001 
Bazell & Dwek |1990[). 



The Galactic interstellar extinction curve displays a 
"2175 A peak" , which has presented an astrophysical puz- 
zle since its discovery by Stecher ( |1965[ ). Fitzpatrick & 
Massa ( |1986D studied the interstellar extinction in the di- 
rection of 45 reddened stars, and found that it displays 
a peak whose central wavelength Ao is remarkably con- 
stant (Ao = 2174.4 ± 17 A), even though its full width at 
half maximum (FWHM) varies considerably from 360 to 
600 A; they also found no apparent correlation between 
the small variation in Aq and the large variation in the 
FWHM. 

Graphite has long been considered a very promis- 
ing though controversial candidate f or ex plaining the 
2175 A peak (e.g. St echer & Donn 1965 ; Fitz patrick 



fc M assa |l98^ ; Mathis |l994t V oshchinnikov |l99C|; SorcU 



199C: Draine & Malhotra 1993; Rouleau et al. 1997; Will 



& Aannestad |1999|) . It is considered a promising candi- 
date because small, 0.015 //m, graphitic spheres display 
an extinction peak at about the right wavelength with a 
FWHM in accordance with the observations (Gilra 1972| ) , 
and because its abundance do not seem to contradict the 
cosmic abundance constraints (Snow & Witt 1995| ). It 
is considered a controversial candidate, however, because 
increasing the size of a small graphitic particle simulta- 
neously increases the central wavelength and FWHM of 
its extinction peak. The peak will shift to longer wave- 
lengths when the particle size is increased, when the par- 
ticle shape is oblate spheroidal and when the particles are 
coated with a dielectric substance such as ice; it will shift 
to shorter wavelengths when the particle shape is prolate 



spheroidal (Gilra 1972; Hecht 1981; Draine & Malhotra 



1993). The observed lack of correlation between the posi- 



tion and FWHM of the peak therefore presents a challenge 
to the hypothesis that graphite particles originate this 
peak. In an extensive investigation, Draine & Malhotra 
(|1993D conclude that if graphite particles are the carriers 
of the 2175 A peak, a variation in their optical proper- 
ties ought to be present. This may be the result of vary- 
ing amounts of impurities, variations in crystallinity, or 



changes in its electronic structure due to surface effects. 
Clustering effects have also been considered; Rouleau et 
al. (1997) have shown that compact clusters of graphitic 
spheres satisfy the criterion that the position of the peak 
remain stable while the width varies. 

To investigate the clustering effect further we here 
compute and analyse the extinction of different polycrys- 
talline graphitic clusters. We have chosen clusters ranging 
from small to large, and that are either sparse or com- 
pact. In this way we intend to evaluate how the extinction 
is influenced by the structure. We focus on clusters con- 
sisting of 4, 7, 8, 27, 32, 49, 108 and 343 touching poly- 
crystalline graphitic spheres. The extinction of the clus- 
ters is calculated by the use of a rigorous (GA) method 
(Gerardy & Ausloos |l982|) as we ll as by two DDA meth- 
ods - MarCODES (Markel |l998|) and DDSCAT (Draine 
& Flatau 1994 ) - to test how well these approximations 
perform when applied to clusters with different morphol- 
ogy. Recently Xu & Gustafson (1999) compared light- 
scattering calculations by a rigorous solution similar to the 
GA method, with DDSCAT for two identical polystyrene 
spheres in contact. They found that DDSCAT worked rea- 
sonably well on small volume structures while its validity 
is challenged on large structures. We find a tendency for 
both the DDA codes to overestimate significantly the ex- 
tinction for small (N < 10) fractal clusters and to a much 
larger extent than for the comparable compact clusters. 

By studying clusters consisting of polycrystalline 
graphitic particles we deal with the anisotropy of graphite 
in a different way than usual (e.g. Draine 1988) and it 
turns out that this has an influence on the comparison 
of the calculated extinction with the observed interstellar 
extinction. 

This paper is structured as follows. Section || deals 
with the anisotropy of graphite. Section |^ introduces the 
fractal and compact clusters studied in this work. Section 
^ describes the computational methods under scrutiny; the 
exact GA solution (Gerardy & Ausloos 1982 ) and the two 
different D DA m ethods, developed by Dr aine fc Flatau 
(DDSCAT, |1994| ) and Markel (MarCoDES, |1998[ ). Section 
g| presents and dicusses our results for clusters from one up 
to 343 particles, in the wavelength range 0.03 — 100 ^m. 
Section discusses the implications of our results in the 
interpretation of the interstellar extinction curve. Section 
presents the conclusions. 

2. Graphite 

Theoretical computation of absorption and scattering by 
graphitic particles is difficult because graphite is a semi- 
metal with high anisotropy. Graphite can be characterised 
by two different dielectric functions, £± and e|| , corre- 
sponding to the electric-field vector E being perpendicular 
(ej^) and parallel (e||) to the symmetry axis of the crystal 
(c-axis), which is perpendicular to the basal plane. It is a 
lot easier to experimentally determine e± than e|| , because 
graphite cleaves readily along the basal plane and hence 
reflectivity measurements can be made with normally in- 
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Table 1. The clusters presented in this paper have three 
different geometries: fractal (frac; D= 1-77), face-center 
cubic (fee) and simple cubic (sc). 
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Fig. 1. The average dielectric function {^ave — "s^-L ~^ s^ll' 
see text) of graphite based on the dielectric functions de- 
rived by Draine & Lee (1984). 



cident light, whereas it is very difficult to prepare suitable 
optical surfaces parallel to the c-axis. 

In this work we deal with the anisotropy of graphite by 
assuming that in all our clusters, each individual particle 
is polycrystalline having a dielectric function eavc given 
by the arithmetic average of e|| and e^, namely tavc = 
ie|l + |e^. For a polycrystal, this arithmetic average is in 
fact an attainable upper bound for its dielectric function 
(Avellaneda et al. 1988 ). We use the d ielectric functions e|| 
and e_L derived by Draine & Lee ( 1984 ) covering the region 
from the far-IR to the far-UV. The data agree well in the 
0.03—62 /i m reg ion with those published by Borghesi & 
Guizzetti ( |1991| ). 

In constrast, the usual "1/3—2/3" approximation, 
treats individual particles as mono-crystalline - 1/3 of the 
cluster particles are assumed to have dielectric function 
e|| and the remaining 2/3 to have dielectric function e^. 
This approximation has been shown by Draine (1988) and 
Draine & Malhotra ( 1993| ) to have a surprisingly good ac- 
curacy for graphite grains with radii < 200 A. However, 
considering grain formation and grain growth in stellar en- 
vironments (Sedlmayr 1994), mono-crystalline particles do 
not seem to be as valid an assumption as polycrystalline 
ones. Also, investigations of presolar nano-diamond grains 
from meteorites - direct specimens of surviving physical 
material formed in past stellar environments which can be 
quantitatively analysed in the laboratory - show a clear 
tendency towards being polycrystalline (Daulton et al. 



1994; Phelps 1999) 



A particle in vacuum will show a Lorenz-Mie absorp- 
tion peak whenever the real part of its dielectric function, 
3?(e), satisfies 3?(e) = —2. Looking at the average dielec- 
tric function of graphite shown in Fig. ^ absorption peaks 
should occur at about 0.220 /im and 14.6 /xm. The latter 
should be much more damped due to the higher value of 
the imaginary part of e. 



3. Structure of the clusters 

We consider three-dimensional clusters of identical touch- 
ing spherical particles arranged in three different geome- 
tries: fractal, simple cubic, and face-centred cubic. The 
structures do not have shapes expected to be found in 
space, but will provide us with boundary conditions for 
the problem of calculating the extinction from clusters of 
grains of different morphology. Tablejl] lists all the clusters 
used in this work. 

The fractal clusters, frac7, frac49, and frac343, are ob- 
tained from the first three stages of the recursive con- 
struction of the snowflake fractal. This construction can 
be summarised as follows: a seed particle is put at the 
origin of the coordinate system. In the first step a genera- 
tor is built by symmetrically gluing to the seed particle six 
copies of itself along the x, y, and z axes. Next, each parti- 
cle in the first configuration, the generator, is substituted 
by the whole generator itself. In the next steps the same 
rule is applied: each particle is replaced by the generator. 
Fig. H shows this procedure up to the second stage. For a 
fractal cluster its dimension D is defined by 



M(r) = Mo 



(1) 



where M(r) is the mass of material contained within 
a sphere of radius r. For a solid grain of constant 
density Z? = 3, whereas for a fractal grain _D < 3 
(Mandelbrot 1983| ) ; in parti cular, for the snowflake fractal 
D = In7/ln3 = 1.77 (Vicsek |l98Ej) . 

Although such a deterministic structure (Fig.||) is not 
expected to occur in nature, its fractal dimension is close 
to that of more realistic random cluster-cluster aggrega- 
tion models. In particular, small particles in space may 
move in straight, ballistic trajectories and form larger ag- 
gregates upon collisions. Numerical simulation of this pro- 
cess yields a value of around 1.9 for th e frac tal dimen- 
sion of the r esulting aggregates (Mea kin 198g ; Botet and 



JuUicn |1988| Meakin & JuUien |1988[ ); this value has also 
been confirmed by experimental results (Wurm & Blum 
1998 ). Moreover, optical properties of fractal clusters with 
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Fig. 2. Stages of fractal construction. The seed particle, 
the generator (quasi-fractal clusters consisting of 7 parti- 
cles) and the 3-D structure of the 49 particle fractal cluster 
(D=: 1.77), see text for more details. 



D < 2 are predicted to be significantly different from those 
with D > 2 (Berry & Percival 1986); hence we consider 
important to use a fractal with a realistic dimension in 
our computations. The open fractal structure presents a 
challenge for modelling due to the high porosity and large 
surface area. 

As a contrast to the fractal structure, compact crys- 
talline structures are studied, namely, face-centred cubic 



and simple cubic, see e.g. Kittel (1986) for a discussion 
on crystal structures. All the clusters listed in Table|l| are 
symmetric except the fcc49[|. The non-fractal clusters are 
a challenge to model since they are far from being spher- 
ical. 



4. The computational methods 

4.1. Gerardy and Ausloos theory 

A rigorous and complete solution to the multi-sphere light 
scattering problem has been given by Gerardy & Ausloos 
(GA) ( 119801 ; |l982t |1983|; 



as an extension of the Mie- 



Ruppin theory (Mie 1908 ; Ruppin 1975 ) . It is based on the 
exact solution of Maxwell's equations for arbitrary clus- 
ter geometries, polarisation and incidence direction of the 
light. This is done by expanding the various fields involved 
in terms of the vector spherical harmonics (VSH). The 
usual boundary conditions are extended to take into ac- 



^ For such an fee cluster 48 particles would form a 3 x 2 x 2 
unit-cell, while 49 particles can only be fitted into a 3 x 2 x 3 
unit-cell. 



count the possible existence of longitudinal plasmonsg in 
the spheres. High-order multipolar electric and magnetic 
interaction effects are included. We summarise their work 
as used throughout the paper. 

We consider a cluster of N homogeneous spheres of 
radius R and dielectric function e(w), embedded in a ma- 
trix of dielectric constant em and submitted to a plane 
polarised time harmonic electromagnetic field. The total 
scattered field from the cluster is represented as a super- 
position of individual fields scattered from each sphere. 
The electromagnetic field impinging on each sphere con- 
sists of the external incident wave and the waves scattered 
by the other spheres. For any sphere, the incident, inter- 
nal and scattered fields are expressed in VSH centred at 
the sphere origin. The boundary conditions on its surface 
are solved by transforming all relevant field expansions 
into the sphere coordinate system, yielding the following 
system of linear coupled equations for the expansion coeffi- 
cients and df^ of the scattered field (Gerardy & Ausloos 



19821) 



df, ^ Ag{ bo ^ + + Jf^i' diy)}, (3) 

where /i — {q,p,i) and v = {n,m,j). The indices q and 
n denote the polar order, p — ~q, . . . ,q, m = — n, . . . , n, 
and the indices i,j number the particles in the cluster. 
Moreover, oq^ and h^^ are the coefficients of the expan- 
sion of the external incident wave in VSH in the coordi- 
nate frame of the i th sphere (Gerardy & Ausloos 1982). 



The multi-polar electric and magnetic susceptibilities of a 
sphere (Gerardy & Ausloos 1982| ) are denoted by Ag and 
Tq respectively. The coefficients J^^ and C^y describe the 
transformation of the VSH from a frame centred on par- 
ticle j to another centred on particle i. Analytical expres- 
sions for these coefficients have been given by Gerardy and 
Ausloos (|l982|) . The primes in the sums in Eqs. (||) and ^ 
indicate that terms with j = i are omitted. The solution 
to this system of equations can be written in matrix form 
as 



c 
d 



Tjvil TAri2 
TjV21 TAr22 



ao 
bo 



(4) 



where Tjv is the T-mat rix (W aterman 1971; Tsang et al. 
1985 ; Mishchenko et al. 20001 ) of the cluster with compo- 
nents 



(5) 



Tnii = ((r-i-j)-c(A-i-j)-ic)-i 

TAri2 — TjvilC(A ^ — J) ^ 
TAr21 = (A ^ — J) ^CTatii 
Tn22 = (A-l-J)-l(CTAri2+I), 

and I is the identity matrix. Whe n N = 1 , Eq. (|4|) is just 
the Mie- Ruppin result (Mie |l908t Ruppin |l975|) 



c 




r 







ao 


d 







A 




bo 



(6) 



A collective excitation for quantized pleisma oscillations, in 
which the free electrons in a metal are treated as a plasma. 
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The extinction efficiency, Qexti i-e. the extinction 
cross-section in units of the total geometrical cross-section 
can then be obtained from Gerardy & Ausloos (1982) 



Qe 



1 



5R(a5^c + b5'd 



(7) 



where k = 2t:/X and A is the wavelength in the matrix. 
The extinction cross section per unit volume is then com- 
puted as 



Limiting q and n in Eqs. (|^) and (H) to integers less 
than or equal to L, we obtain a system of 2NL{L -f 2) 
equations whose solution, Eq. (^, is the 2^-polar ap- 
proximation to the electromagnetic response of the clus- 
ter; ao, bo, c and d are complex vectors of dimension 
NL{L + 2), while the T^ij, i,j = 1, 2, are complex square 
matrices of dimension NL{L + 2). In this case Eqs. (0) and 

determine the 2^-polar approximation to the cluster 
extinction efficiency and the cluster extinction per unit 
volume respectively. 

The extension of this theory to treat the case of ab- 
sorbing embedding media is given in (Lebedev & Stenzel 
1999| ; Lebedev et al. |l999|) . 



4.2. The DDA method 

The discrete dipole approximation (DDA) - also known 
as the coupled dipole approximation - method is one of 



several disc retisa tion methods (e.g. Draine 1988 ; Hage & 
Greenberg 1990) for solving scattering problems in the 
presence of a target with arbitrary geometry. The dis- 
cretisation of the integral form of Maxwell's equations 
is usually done by the method of moments (Harrington 
1968). Purcell & Pennypacker (1973) were the ffist to ap- 



ply this method to astrophysical problems; since then, the 
DDA method ha s been improved greatly by Draine ( |1988|) , 
Goodman et al. (1991), Draine & Goodman (1993), Draine 
& Flatau ( |1994D , Markcl ( |1998[ ), and Draine ( |2000| ). The 
DDA method has gained popularity among scientists due 
to its clarity in physical principle and the FORTRAN im- 
plementation which have been made publicly available by 
Draine & Flatau (DDSCAT package, Draine & Flatau 
199^ ) and by Markel (MarCoDES, Markel p^98| ). 

When considering the problem of scattering and ab- 
sorption of linearly polarised light 



Eq — cq exp (ikr) 



(9) 



by an isotropic graphitic grain, then within the concept of 
the DDA method, the grain is replaced by a set of discrete 
elements of volume Vi with relative dielectric constant 
and dipole moments di — d(ri), i = 1,...,7V, whose co- 
ordinates are specified by vectors r^. The equations for 
the dipole moments can be written using simple consider- 
ations based on the concept of the exciting field, which is 



equal to the sum of the incident wave and the fields of the 
rest of the dipoles in a given point 



di in) = a. 



Eoin) + k^J2Gv<iAr,) 



(10) 



where the dipole scattering tensor Gy has the following 
form in dyad notations: 



A (kr) % + B (kr) ^ 



(11) 



(8) A{x) ^ [x ^ + ix - X ^] exp {ix) 



B (x) = [-X ^ - Six ^ + 3x^^] exp (ix) 



The solution of the set of Eqs. (|l^) and (0) yields all 
the basic optical characteristics of a particle such as the 
integrated extinction (Qext) and absorption (Qabs): 



Qcxt = 47rfcQ 



/ , (eodj)exp(-ifcorj) 

i 

47r3 (e,) 



(12) 



and scattering (Qsca = Qcxt — Qabs) efficiencies. Here 9 
means the imaginary part of the expression in the argu- 
ment. 

Once the location and polarisability of the points are 
specified, calculations of the scattering and absorption of 
light by the array of polarisable points can be carried out 
to in principle whatever accuracy is required. The limiting 
factor is the capacity of computing resources. 

4.2.1. DDSCAT 

In this work we use the DDSCAT code version 5alO 
(Draine & Flatau |1994| ; Draine & Flatau |2000| ). This ver- 



sion contains a new shape option where a target can be 
defined as the union of the volumes of an arbitrary num- 
ber of spheres. In DDSCAT the considered grain/cluster 
is replaced by a cubic array of point dipoles. The cubic ar- 
ray has numerical advantages because the conjugate gra- 
dient method can be efficiently applied to solve the matrix 
eq uation describing the dipole interactions (Goodman et 
al. 1991). When knowing the dipole strength of each cell in 
the particle it is possible to compute the optical properties 
of arbitrary dust configurations. 

There are three criteria for validity of DDSCAT: 

1. The wave phase shift p = Imlkd (m = ^/e being the 
complex refractive index of the target material) over 
the distance d between neighbouring dipoles should be 
less than 1 for calculations of total cross sections and 
less than 0.5 for phase function calculations. 

2. d must be small enough to describe the object shape 
satisfactorily. 

3. The refractive index m must fulfill ItoI < 2. 
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A comparison study by Draine & Flatau (1994) shows 
that scattering and absorption cross sections can be cal- 
culated with DDSCAT to accuracies of a few percent, pro- 
vided that the criteria elaborated above are satisfied. 

Xing & Hanner ( 1997 ) state that the Q efficiency fac- 



tor is generally much less sensitive to those criteria, and 
that it is usually sufficient to let the value of p be around 1 
or possible bigger. Xu & Gustafson ( 1999|) show through a 
comparison between the exact solution for two spheres in 
contact an d DD SCAT that the criterion set up by Draine 
& Flatau (1994) is not sufficient to ensure high accuracy 
if the particles are large and strongly interacting (high 
refractive index). They, therefore, recommend to use a 
validity criteria of p < 0.3 when calculating total cross 
sections. 

The typical number of dipoles needed to obtain a reli- 
able computational result using the DDSCAT code can be 
determined by calculating the minimum number of dipoles 
needed per particle. When a particle is represented by a 3- 
dimensional array of N dipoles, its volume is Nd^ , which 
must be equal to 47r_R'^/3, 



_ 47r /i?y _ 47r /27ri?|m| 



1039 



R\m\ 



(13) 



since d is related to the wave phase shift phy d — p/{\m\k) 
(Draine & Flatau |200C|; Xing & Hanner [l997|) . The number 



of dipoles needed to satisfactorily describe the shape of 
the object considered might be larger than indicated by 
Eq. ([1^) if the shape is "challenging" . It is important to 
remember that Eq. (|l^) only ful fills t he first of the three 
criteria set by Draine & Flatau (200C). 

For materials with large refractive index (|m| > 2), 
Draine & Goodman ( 1993| ) show that especially the ab- 
sorption is overestimated by DDA. The limitation in 
DDSCAT is set by the use of the lattice dispersion relation 
(LDR) for electromagnetic waves propagating on an infi- 
nite cubic lattice of point dipoles of polarisability and 
spacing d. According to Draine & Goodman ( 1993] ), the 
LDR prescription for ai gives fair accuracy for scattering 
but poorer results for absorption. When |m— 1| is large the 
continuum material is effective at screening the external 
field: in the limit |to — 1| — *■ oo the internal field gener- 
ated by the polarisation would exactly cancel the incident 
field, so that the continuum material in the interior of the 
target would be subjected to zero field. In the case of a 
discrete dipole array, the dipoles in the interior will also be 
effectively shielded, while the dipoles located on the tar- 
get surface are not fully shielded and, as a result, absorb 
energy from the external field at an excessive rate. This 
error can be reduced to any desired level by increasing 
the number N of dipoles, thereby minimising the fraction 
N~^/^ of the dipoles which are at surface sites, but very 
large values of N are required when |m| is large (Draine 
& Goodman [199^ ). 

Graphite is characterised by having a high refractive 
index. This means that the criterion \m\ < 2 is only ful- 
filled shortward of 0.072 fim {m = n + ik = 0.58 + il.42) 



and relaxing the criterion a little (to account for the vari- 
ability of TO that occurs) gives an upper limit of 0.216 /xm 
(to = 0.66 + il.35). 

4.2.2. MarCoDES 

Another efficient code based on DDA is the Markel 
Coupled Dipole Equation Solver (MarCoDES, Markel 



1998 ). This code is designed to solve the coupled Eq. (10) 



for an arbitrary cluster of point dipoles using either con- 
jugate gradient method (iterative method) or the LU ex- 
pansion (direct method). The program is in principle ap- 
plicable to clusters of small particles with arbitrary ge- 
ometry, but is most computationally efficient for sparse 
clusters (i.e. when the volume fraction is very low) with 
significant number of particles (« lO'^ — lO'*). Contrary 
to DDSCAT the program does not use the Fast Fourier 
Transformation (EFT) because this might significantly de- 
crease the computational performance for clusters with a 
low volume filling factor. When the volume filling fraction 
is close to unity, algorithms utilising EFT will be much 
faster. In the framework of the MarCoDES a spherical 
grain is considered to be equivalent to the single elemen- 
tary dipole, this corresponds to iV = 1 in the DDSCAT 
case. The dimensionality of the coordinates of particles 
in MarCoDES require a special consideration. By replac- 
ing real particles by point dipoles located at their cen- 
tres we significantly underestimate the strength of their 
interaction. In order to correct the interaction strength, 
the author of MarCoDES introduces geometrical inter- 
section of particles. All coordinates are defined in terms 
of the distance between neighbouring dipoles d, which is 

1/3 

given hy d = (47r/3) R. So, for example, if two particles 
have radii 10 nm, then the distance between the dipoles is 
16.12 nm. This suggested phenomenological procedure al- 
lows MarCoDES to be more accurate than the usual single 
dipole approximation since the intersection produces some 
analogy of including higher multipole interactions between 
particles. Application of MarCoDES to fractal sparse clus- 
ters allows, in principle, to introduce renormalisation on 
the sphere radii and number of particles depending on the 
fractal dimension of the cluster (Markel et al. 200C), which 
will produce even further corrections to the strength. In 
this paper we use the simple particle intersection model 
given in the MarCoDES. The fact that the program only 
uses a single dipole for each particle in the cluster have sig- 
nificant benefits in computation efficiency when compared 
to other multi-polar approaches such as the GA method 
and DDSCAT. 

5. Results 

We present the results of the computation of the extinc- 
tion of infrared, visible, and ultraviolet light by clusters of 
graphitic spheres of radii 10 nm and 50 nm, containing as 
few as 1 sphere to as many as 343 spheres. Specifically, our 
results concern the wavelength region from 0.1 to 100 /im. 
For convenience, we shall, from now onwards, call a clus- 
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ter small if it has less than 10 particles, medium-size if it 
has between 25 and 50 particles, and large if it has more 
than 100 particles. With this convention, we shall first 
discuss the convergence of the GA computations, followed 
by a comparison of the three computational methods (GA, 
MarCoDES, and DDSCAT) when appUed to the calcula- 
tion of the extinction of a single polycrystalline graphitic 
sphere. Then we shall present in succession our results for 
small, medium-size, and large clusters, and end the section 
with an overall assessment of the computational methods 
in light of the results presented. 

5.1. Convergence of GA method 



As we stated in Sect. 4.1, with the GA method we can 
compute the extinction of a cluster in the 2^-polar ap- 
proximation. We find that the smallest L needed for the 
convergence of the extinction of our clusters will be differ- 
ent for different regions of the optical spectrum. For small 
clusters, our full converged results show that in the UV- 
visible it is sufficient to use L = 7 for open clusters, and L 
= 9 for compact ones, to compute the extinction with an 
accuracy of 1%; whereas for longer wavelengths, at least L 
= 11 is required to ensure the same accuracy; this is so be- 
cause in this region graphite has a metallic-like behaviour 
(high |m|), and, thus, higher order multipoles have to be 
included in the computation of the extinction to get the 
same accuracy as in the UV- visible range. We note that 
in the UV-visible range, by accepting an accuracy of 5% 
in the computation of the extinction, we can use L = 5 for 
open clusters and L = 7 for compact ones; we expect this 
to hold for clusters of up to a few tens of particles. 

In column 2 of Table ||, we give the cut-off polar order L 
used in the calculation of the extinction of all our clusters. 
Full convergence was achieved only for small clusters; for 
the others, L indicates the maximum polar order obtained 
with our computer resources. In general, L = 11 was used 
for small clusters, L from 6 to 7 for medium-size clusters, 
and L from 2 to 3 for large clusters. 

We show in Fig. || the convergence of the extinction of 
the frac7 cluster in the wavelength range 0.03 — 100 fim. 
The inset shows the 2175 A peak; around the peak maxi- 
mum L = 3 already gives results accurate within 1%. 
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Fig. 3. Solution of the GA method for the fractal cluster 
structure containing 7 polycrystalline graphitic particles. 
The solution is shown for different polar orders L. At L = 
11 the solution was fully converged over the whole wave- 
length interval. The particle radius was 10 nm. 
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Fig. 4. A comparison between of DDSCAT, MarCoDES 
and the exact GA solution for a polycrystalline graphitic 
sphere with a 10 nm radius. The GA solution was fully 
converged at L = 2. The DDSCAT calculation was done 
with 24464 dipoles. 



5.2. Single sphere 

To set up a comparison baseline, we compute the extinc- 
tion of af single sphere of radius 10 nm using the two DDA 
codes (DDSCAT and MarCoDES) and the GA solution. 
Since the GA theory extends that of Mie (Sect. 4.1), the 
GA and Mie solution coincide for a single sphere; this so- 
lution will be taken as a reference for the comparison. 

It can be seen from Fig. I that DDSCAT is almost 
indistinguishable from the Mie solution for A < 1.0 /im 
while MarCoDES is so for A > 0.6 /im. The less good fit 
by DDSCAT for A > 1.0 /xm is most likely a consequence of 
the refractive index of graphite having a modulus larger 



than 2 for these wavelengths, and of the 24 464 dipoles 
being inadequate to account for this. 

5.3. Small (N < W) clusters 

As before, we take the fully converged GA solutions as 
a reference for comparison. Fig. || shows the extinction of 
the frac7 cluster as computed with the three methods, 
while Fig.^ shows similar information for the sc8 clus- 
ter. For the frac7 cluster (Fig.||) MarCoDES underesti- 
mates the extinction for 0.2 /im < A < 0.25 fim and sys- 
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Fig. 5. A comparison of the solutions from GA, 
MarCoDES and DDSCAT of the extinction cross section 
of the fractal cluster, shown in Fig. ^ consisting of 7 poly- 
crystalline graphitic particles. The GA solution was fully 
converged at L = 11. In the DDSCAT calculation 6 838 
dipoles were used. The particle radius was 10 nm. 



Fig. 6. A comparison of the solutions from the GA, 
MarCoDES and DDSCAT calculations for the extinction 
cross section of the simple cubic cluster containing 8 poly- 
crystalline graphitic particles. The GA solution was fully 
converged at L = 11. In the DDSCAT calculation 16 824 
dipoles were used. The particle radius was 10 nm. 



tematically overestimates it for A > 0.26 /im. DDSCAT0 
slightly underestimates the extinction for A < 0.23 /im and 
overestimates it for A > 0.23 nm. In the range 0.26 /im 
< A < 4 /im the excess extinction from DDSCAT is 
smaller than for the MarCoDES calculations while for 
A > 4 /tm DDSCAT gives a systematic relative increase 
in the excess extinction with wavelength which is much 
larger than the almost constant factor, ~ 1.5, seen for 
MarCoDES. 

For the sc8 cluster (Fig.|^) the trend is different. 
DDSCAT^ overestimates the extinction for A < 0.3 /im, 
underestimates it for 0.3 /tm < A < 11 /im and overesti- 
mates it again for A > 11 /im. MarCoDES also overesti- 
mates the extinction for A < 0.26 /tm but then system- 
atically underestimates it to a much larger extent than 
DDSCAT for A > 0.26 /tm. 

We have also investigated the effect of particle-size on 
a cluster's extinction using the fully converged results of 
the GA solution. For this, in addition to the extinction of 
the frac7 cluster of spheres of radii 10 nm, we have also 
computed the extinction of a fracT cluster of spheres of 
radii 50 nm; which was fully converged at L = 11 with 
a precision of 1.2%. Fig. ^ shows the results; what strikes 
the eye immediately is the high increase in the extinction 
over almost the entire wavelength range investigated, upon 
increasing the radii of the spheres. Still, for wavelengths 
A > 1.1 /tm the shape of the extinction of both clusters 
look roughly the same. More important, however, are the 
differences observed in the wavelength range A < 1.1 /i, 

^ Calculated with a 36 x 36 x 36 dipole grid, which provides 
6838 dipoles (~ 977 dipoles per particle). 

* Calculated with a 32 x 32 x 32 dipole grid, which provides 
16 824 dipoles (~ 2103 dipoles per particle). 
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Fig. 7. The extinction cross section obtained with GA for 
frac7 with particle radius of 10 and 50 nm. The calcula- 
tions were fully converged at L = 11. 



where the extinction of the cluster with the larger spheres 
exhibits peaks at 0.11 /im and 0.4 /im in addition to the 
2175 A peak, which is the only peak displayed by the 
extinction of the cluster with the smaller spheres. All these 
observations agree with the finding of Kimura (2001) that 
the light scattering of fractal clusters depends strongly on 
the size parameter of the monomer. 

5.4. Medium-size (25 < N < 50) clusters 

For the frac49 cluster, shown in Fig.||, the extinction 
cross section obtained by the three methods is dis- 
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0.17 0.19 0.21 0.24 0.27 0.30 



0.1 1.0 10.0 100.0 

Wavelength (^m) 

Fig. 8. A comparison of the solutions from GA, 
MarCoDES and DDSCAT for the extinction cross section 
of the frac49 cluster, shown in Fig.||. The GA solution 
was truncated at L = 6. In the DDSCAT calculation 1 722 
dipoles were used. The particle radius was 10 nm. 
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Fig. 10. A comparison of the solutions from the GA and 
DDSCAT calculations for the extinction cross section of 
the fcc32 cluster. The GA solution was truncated at L = 
7. In the DDSCAT calculation 17 942 dipoles were used. 
The particle radius was 10 nm. 



- GA, L=7 
MarCoDES 
DDSCAT, N=16783 
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Fig. 9. A comparison of the solutions from the GA, 
MarCoDES and DDSCAT calculations for the extinction 
cross section of the sc27 cluster. The GA solution was 
truncated at L = 7. In the DDSCAT calculation 16 783 
dipoles were used. The particle radius was 10 nm. 



played in Fig.||. For this structure neither DDSCAl]^ nor 
MarCoDES gets close to the solution suggested from the 
GA calculations with L = 6. The GA solution is not fully 
converged at long wavelengths due to computational lim- 
itations; however, in the UV-visible region convergence is 
ensured at least with an accuracy of 5%. Both MarCoDES 
and DDSCAT greatly overestimates the extinction cross 
section for this cluster at longer wavelengths. 



As a comparison to the frac49 cluster we have car- 
ried out calculations for two symmetric clusters, sc27 and 
fcc32, and an asymmetric one, fcc49. 

For the sc27 cluster (Fig. ^ the GA solution was trun- 
cated at L = 7. MarCoDES underestimates the extinction 
for 0.17 fim< X< 0.22 fj.m and A > 0.29 /zm. DDSCAT| 
also slightly underestimates the extinction for 0.17 fj,m 
< X < 0.22 iim and 0.26 /im < A < 0.4 /^m, but system- 
atically overestimates it for A > 1 /im. 

We found that the version of MarCoDES tested here 
can not calculate a face-center cubic structure of touching 
particles because in this case the lattice cells represent- 
ing neighbouring particles will touch only at the corners, 
giving as a consequence a spectrum corresponding to non- 
touching particles. The extinction of fcc32 and fcc49 clus- 
ters were therefore only calculated with the GA theory and 
DDSCAT. The results for fcc32 are displayed in Fig. ^ 
the results for the fcc49 cluster look similar. The GA cal- 
culations were truncated at L = 7 for fcc32 and at L = 
6 for fcc49. For the DDSCAT calculations a 32 x 32 x 32 
dipole grid was used; for the fcc32 cluster this provides 
17942 dipoles (~ 561 dipoles per particle). For fcc49 a 
48 X 32 X 40 dipole grid was used; this provides 29 849 
dipoles (~ 609 dipoles per particle). As can be seen in 
Fig|l^, DDSCAT systematically overestimates the extinc- 
tion especially at long wavelengths. The overestimation of 
the extinction cross section as calculated by DDSCAT is 
a general tendency for all the clusters we have studied. 
Also DDSCAT does not reproduce the absorption hump 
expected for graphite around 10—15 fim, which is repro- 
duced in the GA and MarCoDES calculations and also in 



Lorenz-Mie calculations by Draine & Li (2001) 



^ Calculated with a 36 x 36 x 36 dipole grid, which provides 
1722 dipoles (~ 35 dipoles per particle). 



Calculated with a 32 x 32 x 32 dipole grid, this provides 
16 783 dipoles (^ 622 dipoles per particle). 
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Fig. 11. A comparison of the extinction cross section ob- 
tained with GA theory for the clusters frac49, fcc49, fcc32 
and sc27. The 49 particle clusters were truncated at L = 
6, while the two other clusters were truncated at L = 7. 
The particle radius was 10 nm. 



5.5. Large (N > 100) clusters 

The GA solution could only be calculated to L = 3 for 
the fcclOS cluster and to L = 2 for the frac343 cluster 
due to computational limitations. The extinction of the 
fcclOS and frac343 clusters looks very similar to that of 
the fcc32 and frac49 clusters shown in Figs. ^ and ||, re- 
spectively. For the frac343 cluster the whole spectrum is 
shifted to longer wavelengths by about 0.002 ^m and for 
the fcclOS cluster the whole spectrum is shifted 0.01 /im 
to longer wavelengths compared to the frac49 and fcc32 
cluster, respectively. It was not possible to calculate these 
large clusters with DDSGAT. 



interesting to note that MarCoDES gives better results 
for compact than for fractal clutcrs of intermediate size. 
According to Markel et al. ( ^000 ), the performance of 
MarCoDES can be improved by altering the intersection 
parameter which determines if the particles touch or over- 
lap, but we have not investigated that here. Nevertheless, 
the fact that MarCoDES in some cases overestimates the 
extinction, as compared to GA theory, and in other cases 
underestimates it, suggests that the determination of the 
rather arbitrary optimal intersection parameter is a very 
complex problem indeed. 

We next compare the extinction of fractal and compact 
clusters of intermediate size. In Fig.|ll|, the GA calcula- 
tion for frac49 is compared to that for fcc49 and fcc32. It 
is seen that the extinction of fractal and compact clusters 
are of the same order of magnitude. This is in contradic- 
tion to the DDA calculations which give a much higher 
extinction for fractal clusters (Figs.|5| - ^l|). We, there- 
fore, recommend that great care should be taken when 
using these DDA codes at longer wavelengths for fractal 
clusters with low fractal dimension consisting of materials 
with high refractive index (m). 

Concerning the effect of particle size on the extinc- 
tion cross section, in classical electromagnetic theory the 
extinction is independent of cluster size in the long wave- 
length limit, i.e. when the size is much smaller than the 
wavelength. As seen in Fig.^, clusters consisting of par- 
ticles of radii 50 nm are clearly outside of this limit. 
Decreasing the particle radius below 10 nm, however, will 
give, in a classical theory, at most minor effects on the 
extinction of single particles and small clusters. Quantum 
mechanical effects, on the other hand, may influence the 
width and the position of the absorption peak for such 
small particles. A detailed study of these particle size ef- 
fects on the extinction cross section is outside the scope 
of the present work. 



5.6. Discussion 

It is seen that both DDA codes overestimate the extinc- 
tion of the fractal structures at long wavelengths. A likely 
explanation is that the DDA model considers only point 
dipole interactions, whereas the GA solution matches the 
boundary conditions on the surface of the spherical parti- 
cle and solves the problem of scattered and internal fields 
at the same time. Also, the large surface area of frac- 
tal clusters contributes to the complexity of the scattered 
field, which the GA method handles well because it takes 
into account the surface topology through the boundary 
conditions; the DDA codes, on the other hand, perform 
best when the surface to volume ratio is low. This in- 
terpretation is also supported by the calculations for the 
small compact sc8 cluster, see Fig.||, where the DDA ap- 
proximations perform reasonably well. For larger compact 
clusters, however, DDSGAT does not give good results at 
long wavelengths. This is due to computational limita- 
tions; the number of dipoles per particle that we could 
handle becomes lower as the cluster size increases. It is 



6. The 2175 A extinction peak 

The main observational constraints concerning the 2175 A 
peak are according to Fitzpatrick & Massa (1986) and 
Mathis (|1994|): 



1. The peak position is remarkably constant, fmax = 
4.6 ± 0.04 /im"\ where v = l/X. 

2. The peak width has a much larger range of variations 
of 1.0 ± 0.25 Aim-i (i.e. < 25%). 

3. The variation in peak position and the width (7) 
are uncorrelated, except for the widest bumps (7 > 
1.2 /im~^) where a systematic shift to larger wavenum- 
ber is observed (Cardelli & Savage 1988; Cardelli & 
Clayton 1991). These later observations all have lines 
of sight passing through dark, dense regions. 

Table || lists the peak position and width for the differ- 
ent cluster structures, as computed with the GA method. 
The polar order at which the GA calculation was trun- 
cated is also indicated. The width was determined as the 
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Fig. 12. Extinction cross section for polycrystalline 
graphitic clusters of different sizes calculated with the GA 
method. A difference in peak position can be seen be- 
tween the fractal clusters and compact clusters. The lines 
in the plot indicate the peak position and width variation 
of the interstellar extinction curve. All clusters consist of 
particles with a 10 nm radius. Also shown is the GA solu- 
tion for a sphere, radius 10 nm, which is equivalent to the 
Lorenz-Mie solution. 

Table 2. Peak position and FWHM of different clusters, 
as calculated with the GA method. The value of L indi- 
cates at which polar order the GA calculations were trun- 
cated. 



cluster name 


L 


peak 


FWHM 










frac7 


11 


4.52 


1.17 


frac49 


6 


4.52 


1.33 


frac343 


2 


4.42 


1.27 


sphere 


3 


4.62 


0.96 


fcc4 


11 


4.43 


1.15 


fcc32 


7 


3.98 


1.43 


fee 49 


6 


3.71 


1.52 


fcclOS 


3 


3.98 


2.01 


sc8 


11 


4.41 


1.27 


sc27 


7 


4.12 


1.62 



FWHM. For the asymmetric compact cluster, we have in- 
vestigated the importance of the orientation of the cluster 
and it was found that the peak position was shifted less 
than 0.1 /im^^ while the width of the feature was not 
affected. 



It was shown in a study by Rouleau et al. (1997), that 
small compact clusters qualitatively satisfy the observa- 
tional constraints on the 2175 A peak, except that the 
peak position falls at the wrong wavelength. This result is 
not confirmed from our findings, since the compact clus- 
ters we have studied do not have a stable peak position. 
The clusters considered by Rouleau et al. (1997) all had 
peak positions at higher wavenumber than the 2175 A 



peak while the clusters we have considered all have peak 
position at lower wavenumbers. Rouleau et al. (1997) also 
used optical constants from Draine & Lee ( 1984| ), but they 
dealt with the anisotropy of the graphitic particles with 
the usual "1/3— 2/3" approximation, see Sect. ^ which 
assumes the particles to be mono-crystalline. As argued in 
Sect. I assuming the particles to be polycrystalline seems 
much more reasonable. This implies that variations in 
crystallinity is indeed an important factor when discussing 
the 2175 A peak as also suggested by Draine & Malhotra 



( p93D . 

The fractal clusters considered here do not enhance 
the long wavelength wing and do not introduce additional 
structure as seen by Rouleau et al. (1997) for elongated 
and "fluffy" structures. The frac7 and frac49 clusters do 
come close (peak position is off by 0.04 — 0.08 /im~^) to 
the observational constraints. This indicates that small 
{N ~ 5 — 100) fractal clusters ought to be investigated in 
more detail to determine if fractal clusters with low fractal 
dimension have a stable peak position around 4.6 /im^^ 
and produce a variable width depending on the number 
of particles in the cluster. 



7. Conclusions 

We computed the extinction of clusters of polycrys- 
talline graphitic particles in the wavelength range 0.1 to 
100 fim. Computations by the rigorous multipolar theory 
of Gerardy & Ausloos (1982; GA) were compared to two 
formulations of the discrete dipole approximation (DDA) . 

We have compared the extinction of open fractal clus- 
ters and compact clusters. The fractal and compact clus- 
ters display an extinction at long wavelengths, of the same 
order of magnitude as computed with the GA method. It 
seems that the DDA approximations grossly overestimate 
the long-wavelength extinction of small fractal structures. 
If this also holds true for larger fractal clusters is not pos- 
sible for us to say due to computational limitations. The 
DDA codes should, however, be used with caution for this 
type of problem. 

The GA computations were compared to the observed 
interstellar extinction curve. We found that results for 
small and medium-size (less than 50 particles) fractal clus- 
ters are in fair agreement with observational constraints, 
while those of compact sc and fee clusters are not. 

Convergence of the GA computations for small (less 
than 10 particles) clusters was found by including 11 
multipoles. Good approximative results were obtained for 
medium-size (between 25 and 50 particles) clusters. In 
most cases we found a substantial differences between 
the GA theory and the DDA approximations of Draine & 



Flatau ( |1994| DDSCAT) and Markel (|l998f MaxCoDES). 

As predicted by Draine & Goodman ( |1993 ) the ab- 
sorption is significantly overestimated by DDSCAT at 
the wavelengths where the modulus of the refractive in- 
dex, |m|, is large, see Figs.|8[ || and |o[ For |r7i| < 2, 
DDSCAT performs well and better than the MarCoDES. 
MarCoDES, however, is computationally much faster than 
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DDSCAT and the GA calculations. With DDSCAT the 
accuracy is directly dependent on the number of dipoles 
used in the approximation. Hence, by using more dipoles 
than we have used here, a better accuracy could be ob- 
tained, but then the computational cost would increase. 
Which DDA code is best to use will therefore depend on 
the type of problem one wants to address and the accuracy 
needed. 
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